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The role of the potential energy landscape in determining the relaxation dynamics of model clus- 
ters is studied using a master equation. Two types of energy landscape are examined: a single funnel, 
as exemplified by 13-atom Morse clusters, and the double funnel landscape of the 38-atom Lennard- 
Jones cluster. Interwell rate constants are calculated using Rice-Ramsperger-Kassel-Marcus theory 
within the harmonic approximation, but anharmonic model partition functions are also considered. 
Decreasing the range of the potential in the Morse clusters is shown to hinder relaxation towards 
the global minimum, and this effect is related to the concomitant changes in the energy landscape. 
The relaxation modes that emerge from the master equation are interpreted and analysed to ex- 
tract interfunnel rate constants for the Lennard- Jones cluster. Since this system is too large for a 
complete characterization of the energy landscape, the conditions under which the master equation 
can be applied to a limited database are explored. Connections are made to relaxation processes in 
proteins and structural glasses. 
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I. INTRODUCTION 



Some of the most interesting processes in chemical 
physics involve relaxation from a non-equilibrium state. 
Examples include the folding of a protein from a dena- 
tured conformation and the formation of a crystal or glass 
upon cooling a liquid. Given a method for calculating the 
rate constants for processes between mutually accessible 
states, the evolution of a non-equilibrium probability dis- 
tribution can be described by a master equation jlj . 

A natural way to define a state in the master equa- 
tion is provided by an "inherent structure" analysis of 
the potential energy surface (PES) [||. Except at high 
temperatures, the configuration of an interacting system 
oscillates in the basin of attraction surrounding a local 
minimum on the PES, and sporadically undergoes tran- 
sitions into neighbouring basins of attraction. A local 
minimum can therefore be regarded as a single state in 
the master equation, and transition states on the PES 
provide the means for dynamics to occur between the 
minima. We have recently obtained databases of min- 
ima and transition states for a variety of systems , 
providing the necessary ingredients for a master equation 
study. There are at least two advantages to modelling re- 
laxation in this coarse-grained state-to-state way. Firstly, 
the master equation can usually be solved for much longer 
time scales than are accessible by direct simulations in 
which the equations of motion are integrated. Secondly, 
the master equation describes the relaxation of an ensem- 
ble without the need for explicit averaging over separate 
trajectories. In fact, the master equation can be used 
as a guide for simulations, for example to devise optimal 
annealing schedules 0. 



In this contribution, the master equation is applied to 
structural databases that we have previously derived for 
some atomic clusters. Section |fj summarizes the master 
equation technique and methods for obtaining state-to- 
state rate constants in the microcanonical and canoni- 
cal ensembles. Section III presents results for 13- atoms 
Morse clusters, M13. The energy landscapes of these clus- 
ters each resemble a funnel, in which the minima are or- 
ganized into pathways of decreasing energy that lead to 
the global minimum on the PES. The characteristics of 
the funnel change with the range of the potential, and 
have been studied in detail in previous work Sec- 
tion III also includes a discussion of harmonic and an- 



harmonic partition function models for describing equi- 
librium properties of the clusters. In Sec. IV, dynamics 
on the paradigmatic double-funnel energy landscape of 
the 38-atom Lennard- Jones cluster, LJ38, are studied. 
We have previously made a number of predictions con- 
cerning the dynamics of M13 Q and L J38 H , which can 
now be examined. We will also present some new ways 
to interpret solutions of the master equation and extract 
information from them. The conditions under which the 
master equation treatment of interwell dynamics is valid 
will also be addressed, especially in the case of LJ38, 
where the knowledge of the energy landscape is incom- 
plete. Finally, Sec.]y| : 
from this work. 



summarizes the main conclusions 
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A. The Master Equation 

Let P(t) be a vector whose components, Pi(t) (1 < 
i < "-min), are the probabilities of the cluster residing 
in a potential well of the geometrical isomer i at time t, 
the total number of such isomers being n m j n - The time 
evolution of these probabilities is governed by 



dPj(t) 
dt 



(1) 



where fey is the first order rate constant for transitions 
from well j to well i. We can set up a transition matrix 
matrix, W, with components 



m—1 



(2) 



so that the diagonal components Wa contain minus the 
total rate constant out of minimum i. This definition 
allows us to write the set of coupled linear differential 
equations (pj) — the "master equation" — in matrix form: 



dP(t) 
dt 



= WP(i). 



(3) 



If W cannot be decomposed into block form then the 
system has a uniquely defined equilibrium state, P eq , for 
which (dP/dt)|p=peq = 0, i.e., W has a single zero eigen- 
value whose eigenvector is the equilibrium probability 
distribution. W is asymmetric, but can be symmetrized 
using the condition of detailed balance: at equilibrium, 



w^pp = w^pp, 



(4) 



so that W t3 = (Pj q /P° q ) 1/2 W t3 is symmetric. W and 

W have the same eigenvalues, Ai, and their respec- 
tive normalized eigenvectors uW and uW are relate d by 
uS 1 ) = Su^ , where S is the diagonal matrix Su = \J P° q . 
Hence, individual components of the eigenvectors are re- 



lated by u 



m) 



/peq The solutiori of Eq (|) is then 



'PP «i W) e A,t 




(5) 



where a™ is component m of . 

Apart from the zero eigenvalue, all the Xj are negative 
jjj . We label eigenvalues and eigenvectors in order of de- 
creasing algebraic value of the eigenvalue, so that Ai = 0, 
and Xj < for 2 < j < n min . As t — > oo, only the j = 1 
term in Eq. (JsJ) survives, and P(t) — > P cq . This limit 
defines the baseline to which the remaining modes decay 
exponentially The size of the contribution of mode j to 
the evolution of the probability of minimum i depends on 
component i of eigenvector j, and on a weighted overlap 



between the initial probability vector and eigenvector j, 
i.e., the term in square brackets in Eq. (||). The sign of 
the product of these two quantities determines whether 
the mode makes an increasing or decreasing contribution 
with time. Combinations of modes with different signs 
give rise to the possibility of the accumulation and subse- 
quent decay of transient populations as probability flows 
from the initial state to equilibrium via intermediates. 

Eq. (||) requires the diagonalisation of the matrix 
W, whose dimension is the number of minima in the 
database, n m ; n . The time required to compute the eigen- 
vectors scales as the cube of the dimension, and the stor- 
age requirements scale as its square. (Although W may 
be sparse, its eigenvectors are not.) However, once di- 
agonalisation has been achieved, Pi(t) can be calculated 
independently for any minimum i at any instant t. The 
only restriction on t comes from the accuracy to which 
the eigenvalues, Ai, can be obtained; if the error is of the 
order 6X, Eq. (||) may diverge as t approaches 1/SX. 

An alternative way of solving the master equation is 
to integrate Eq. (||) numerically. This approach has the 
advantage of not requiring diagonalisation of W, and is 
therefore the only way to proceed for large databases. 
However, it has a number of disadvantages. Firstly, 
knowledge of the eigenvalues and eigenvectors of W is 
useful in interpreting the time evolution of P(i). Sec- 
ondly, accurate integration over long periods can be very 
slow, since the accumulation of numerical error can cause 
the sum of the probabilities to diverge rapidly. Thirdly, 
the full probability vector P(t) (rather than selected com- 
ponents) must be propagated, and the integration must 
start from the time at which the initial probabilities are 
specified. 

If the initial probability vector, P(0), is strongly non- 
equilibrium, many components of P(f) change rapidly as 
soon as the integration starts, and then relax more slowly 
towards equilibrium. Therefore, when numerically inte- 
grating the master equation, the step size required for 
a given accuracy is usually smaller when t is closer to 
zero, and can be enlarged as t grows. To take advantage 
of this, the numerical integration in the present work 
was performed using a Bulirsch-Stoer algorithm with an 
adaptive step size ||. Results from this method coin- 
cided precisely with those of the analytic solution, where 
the latter could be determined. 

The linearity of the master equation rests on the as- 
sumption that the underlying dynamics are Markovian. 
The probability of the transition i — ► j must not depend 
on the history of reaching minimum i, so that the ele- 
ments of the transition matrix are indeed constants for a 
given temperature or total energy. For this restriction to 
apply, states within a potential well must equilibrate on 
a time scale faster than transitions to different minima, 
so that the transitions are truly stochastic. Previous re- 
sults for other clusters [ ^0|JTl| suggest that intrawell equi- 
libration is quite rapid. The Markovian requirement will 
impose an upper limit to the temperatures at which the 
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master equation can be applied to transitions between 
minima, since at high temperatures the phase point does 
not reside in any one well long enough to establish equi- 
librium within it. Division of configuration space into 
the basins of attraction surrounding the minima is less 
useful in this dynamical regime. 



B. Rate Constants and Equilibrium Properties 

To model the probability flow within a database of 
minima using the master equation requires knowledge 
of the rate constants, fey, for their interconversion. For 
the path from minimum j through a particular transition 
state, denoted f , a general form for the rate constant is 
provided by RRKM theory [p| 



wHe) 

hQj(E)' 



(6) 



Here, flj(E) is the density of states associated with min- 
imum j, and 



W*(E)= ( E n\E')dE' 



(7) 



is the sum of states at the transition state with the re- 
active mode removed; fl^(E) is the density of states at 
the transition state excluding this mode, and is the 
potential energy of the transition state. The total rate 
constant, fey, for the process j — ■+ i, is then obtained by 
summing Eq. (^) over all transition states linking j and 
i. 

Rates in the canonical ensemble, i.e., as a function of 
temperature rather than energy, can be obtained from 
Eq. (Q) by Boltzmann weighting thus: 



k](T) 



^k]{E)n {E)e~ E ^ T dE 



(8) 



where Vj is the potential energy of minimum j. Since the 
Laplace transforms 



and 



/>oo 

Zj (T) = / n j (E)e' E/kBT dE 

poo 

Z\T) = / &(E)e- E/kBT dE 



(9) 



(10) 



are the vibrational partition functions of the minimum 
and transition state respectively, we obtain 



k)(T) 



k B TZ^T) 
h Zj (T) ' 



(11) 



Assigning a density of states or partition function to 
each minimum implies that the equilibrium properties of 



the system are described by the superposition principle 
|13|-[l5[ . The total density of states and partition function 
are given by 

n m ; n n m i n 

n{E) = Mi( E ) and z i T ) = z i{ T )- (12) 

i=l i=l 

Application of detailed balance to Eqs. @ and ( |ll|) sup- 
plies the equilibrium occupation probabilities in the mi- 
crocanonical and canonical ensembles: 

P^(E) = ni(E)/n(E), PP(T) = Zi(T)/Z(T). (13) 

These components emerge from the master equation as 
t — > oo. 

The use of fli(E) and Zi(T) to calculate rates for the 
master equation when P ^ P cq [i.e., when Eqs. dl3| ) 
do not hold] still requires thermal equilibration of phase 
points within each potential well. This condition is ful- 
filled under the time scale separation of in ter- and in- 
trawcll motion mentioned at the end of Sec. II A. 



III. MORSE CLUSTERS 

In this section, we apply the master equation to 
databases derived for the 13-atom Morse cluster, M13. 
The Morse potential |l6| can be written as 

V = % V H = e p( - 1 - r '^ r ^ [gPC 1 -^/^) - 2]e, 

(14) 

where r c and e are the dimer pair separation and well 
depth respectively, p is a dimensionless parameter which 
determines the range of the interparticle forces, with 
small p corresponding to long range. Physically mean- 
ingful values include p — 3.15 and 3.17 for sodium and 
potassium JTtJ up to about 14 for Cgo molecules |l8| , |l9| . 
When p = 6, the Morse potential has the same curvature 
as the Lennard-Jones potential at the minimum, e, r c 
and the atomic mass, m, define the system of reduced 
units for the Morse potential independently of p, giving 
(mr-g/e) 1 / 2 as the reduced unit of time. 

Details of databases obtained for four values of p have 
been given previously Q. They demonstrate that the 
energy landscape of these clusters is a single funnel — 
a collection of kinetic pathways leading to a particular 
low-energy structure. However, the funnel has a steeper 
energetic gradient and lower downhill barriers when the 
range of the potential is long. The number of stationary 
points increases dramatically as the range is decreased, 
and the system must then overcome more barriers to 
reach the global minimum. These effects are expected to 
hinder relaxation towards the global minimum. In Sec. 
A, we examine the performance of some partition func- 
tion models for the equilibrium properties of the clusters, 
before addressing the dynamics in Sec. B. 
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A. Partition Function Models 



(16) 



The ingredients for the equations of Sec. II B are the 



densities of states for each minimum and transition state. 
If the harmonic approximation is applied to the vibra- 
tional density of states of a one-component iV-atom clus- 
ter, the result is 



2M (E-Vi)*- 1 
hf G (k - 1)! QiViY' 



(15) 



where v is the geometric mean frequency, k = 3N — 6 is 
the number of vibrational degrees of freedom and hf G is 
the order of the point group of the isomer. The corre- 
sponding partition function is 



Eqs. (|l5|) and ( |l6[ ) can be applied to transition states by 
excluding the reactive mode, in which case k = 3N — 7. 

The harmonic treatment is approximate in two ways. 
Firstly, the anharmonicity of the potential is ignored. 
Secondly, the density of states around each minimum is 
modelled as the surface area of a hyperellipsoid in phase 
space, and the overlap of the hyperellipsoids belonging to 
different minima is neglected. Both these effects become 
more pronounced as the total energy increases; highly 
anharmonic barrier regions are reached, and the phase 
space hyperellipsoids become larger. 
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FIG. 1. Equilibrium occupation probability of the global minimum (top row) and canonical caloric curves (bottom row) 
for Mi3 with p — 4 (left column) and p = 6 (right column). Circles: canonical MC simulations, solid lines: harmonic 
approximation, dotted lines: MB method, dashed lines: MB(t?p = 0.1) method (see text). 



A simple test of the partition function model is to com- 
pare the equilibrium probability of the global minimum 
predicted by Eqs. (|l|) with the fraction of quenches to 
this structure in the course of a simulation. We have 
observed elsewhere fTT| that the harmonic approxima- 
tion works quite well for the equilibrium probabilities of 
LJ7 isomers in the microcanonical ensemble, even in the 
liquid-like regime. Figures |l|a and |l]b show how the har- 



monic approximation performs for the global minimum 
of the M13 clusters with p = 4 and 6 in the canonical en- 
semble, using the databases obtained previously The 
graphs also show the equilibrium probability of the global 
minimum obtained by quenching from canonical Monte 
Carlo (MC) simulations. In the simulations, a spherical 
container of radius 3 a was imposed to prevent evapora- 
tion. At each temperature 2x 10 6 single-particle warm-up 
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steps were performed before collecting data over 2 x 10 8 
steps, quenching every 10 4 steps to find the local min- 
imum in which the phase point resided. The harmonic 
curves depart from the MC results soon after the clus- 
ter begins to melt. For p = 6 (Fig. |l|b), the harmonic 
curve has the correct qualitative shape, but is shifted to 
higher temperature. For p = 4 (Fig. |l]a), however, the 
harmonic approximation underestimates the MC result 
increasingly badly above about k B T = 0.3 e. 

The modelling of partition functions of individual min- 
ima for use in the master equation has been the subject 
of an extensive study by Ball and Berry |2^,^l|. These 
authors considered a variety of analytic forms for Zi(T) 
which attempt to improve on the shortcomings of the 
harmonic approximation. The greatest improvement was 
produced using an anharmonic model based on a first 
order expansion of the density of states for the Morse 
potential f22^]23]| . The anharmonic correction for each 
minimum was derived from the heights of the barriers 
connected to it. The resulting partition function for this 
"Morse barrier" (MB) method is 



zf v) = zf° n f 1 



T 



(17) 



In this expression, rii is the number of transition states 



connected to minimum i, and a 



(') 



(2AV- 



is the 



anharmonicity parameter derived from the barrier height, 
AVj , of the jth transition state connected to minimum 
i. The power a, = min[l, n/rii] compensates for the pos- 
sibility that there may be more than k transition states 
connected to a given minimum, since the original model 
on which Eq. ( |17| ) is based associated one barrier with 
each normal mode. The MB model experiences problems 
when minima are connected by low barriers because the 
full series for the density of states of a Morse oscillator, 
from which the partition function is derived, diverges as 

k B T approaches AVj. This behaviour correctly corre- 
sponds to dissociation of the Morse oscillator, but for 
a cluster isomerisation the density of states should re- 
main finite above the barrier. Ball and Berry |20| sug- 
gested two ways to circumvent this difficulty. In the first 
approach, anharmonic corrections were applied only to 
low-energy minima, which tend not to have small barri- 
ers. This prescription is not entirely satisfactory because 
the low barriers of high-energy minima are a signature of 
the strong anharmonicity in these states, which is then 
ignored by this method. Furthermore, some criterion for 
selecting minima for anharmonic corrections is required. 
The second method, called MB(?7p), involved limiting 
the anharmonic correction a^k B T for each barrier to a 
plateau value, rjp, at temperatures where it would exceed 
this value, i.e., 



zf B W(T) = Z, 



HO 



n 



r] Pl a^k B T 



(18) 



Although a reasonable value of r/p might be expected 
to be around 0.5 (i.e., where k B T = AVj ), the best 
agreement with quenching from constant temperature 
molecular dynamics (MD) simulations was achieved for 
77P « 0.1. 

The global minimum equilibrium probabilities of the 
M13 clusters with p = 4 and 6 using the MB and MB(?yp) 
formulations are shown by the dotted and dashed lines 
in Figs. [l]a and |l|b. The inadequacy of the unconstrained 
MB model is obvious; the large anharmonic corrections 
of the high-energy minima cause these states to dominate 
as soon as they become energetically accessible, making 
the probability of the global minimum plummet at an 
artificially low temperature. For p = 4, the MB(?7p) 
model produces a marginally better match to the MC 
data at low temperature than the harmonic approxima- 
tion, but then deviates more rapidly. For p = 6, quite 
a reasonable improvement is achieved. However, the po- 
sition of the curve is continuously adjustable from the 
MB method (effectively rjp = oo) to the harmonic ap- 
proximation (jip = 0), so it is difficult to argue that the 
improvement is based on physical insight. 

The effect of the partition function model on the equi- 
librium probabilities of minima other than the global 
minimum is hard to gauge for a system where there 
are so many of them, since the quench statistics are 
poor. However, an impression of the overall description 
of the PES can be gained through thermodynamic prop- 
erties derived from the superposition method. The in- 
ternal energy can be obtained from the standard relation 
U = k B T 2 {d In Z/dT)Ny. The harmonic approximation 
yields the classical equipartition result 



m J2 Zf° (V t + nk B T) 



z ao 



(19) 



and the MB(?yp) model [Eq. @] gives 



(/MB(i) P 



1 



i=l 



MB(t7p) 



Vi + nk B T 



+ ai (k B T) 2 J2 

3 = 1 



l + afk B T 



(20) 



where the step function signifies that the derivative of 
the anharmonic correction from a given mode to the par- 
tition function is zero once the plateau has been reached. 
The derivative is also assumed to vanish at the point 

(i) 

a,j k B T = ?7p, where the plateau starts, even though it 
is really discontinuous. 

Figures ||cand0d compare the canonical caloric curves 
of Mi 3 with p = 4 and 6 given by Eqs. @ and @. Also 
shown are the results from MC simulation, obtained by 
adding the kinetic contribution nk B T/2 to the configu- 
rational part given by the simulations. For both clus- 
ters, the harmonic approximation underestimates the in- 
ternal energy. The MB method incorrectly predicts a 
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sharp transition at the temperatures where the equilib- 
rium probability shifts from the global minimum in Figs. 
|l|a and |b. The MB{r) P = 0.1) method, however, repre- 
sents only a marginal improvement on the harmonic ap- 
proximation. We have also found this to be the case for 
the L Jg cluster studied by Ball and Berry [^o| , using the 
model that gave the greatest improvement for the equi- 
librium probabilities of the various isomers. This result 
shows that optimizing the models for the probabilities 
does not necessarily produce the correct thermodynamic 
behaviour. Doye and Wales found that, whilst first or- 
der anharmonic corrections were enough to improve the 
thermodynamic description of LJ55 in the superposition 
method, second order corrections were necessary for the 
smaller LJ i3 ||. 

Attempts to model the density of states of individual 
minima using an analytic function of the energy implic- 
itly assume that the "shape" of the basin of attraction is 
simple enough to be described adequately by such a form. 
In practice, basins of attraction are probably highly com- 
plex objects. For example, in previous work |4[] we saw 
that increasing the range of the potential removes locally 
stable minima, but remnants of these features are likely 
to persist as shoulders or inflections on the PES, so that 
regions of configuration space that were associated with 
shallow minima for a shorter-ranged potential become 
formally associated with other minima when the range is 
increased. These features might explain why the model 
partition function results in Fig. [l] are worse for p = 4 
than for p = 6. 

To illustrate this effect, we have performed micro- 
canonical MD simulations of M43 with p = 4, and LJ13, 
which closely resembles M13 with p — 6. Periodic quench- 
ing was applied, and the Euclidean distance, D, was cal- 
culated between the configuration point taken at the start 
of the quench and the local minimum to which it con- 
verged. The distribution of D for the subset of quenches 
that led to the global minimum is plotted in Fig. || at 
two energies for each cluster. In each case, the lower 
energy has been chosen just above that at which the tra- 
jectory can escape from the global minimum, so that over 
95% of quenches return to the global minimum, and the 
cluster is exploring a large proportion of the catchment 
basin of this structure. The higher energy was chosen 
such that about half the quenches return to the global 
minimum. The corresponding distributions are broader 
and peak at a higher value of D, as would be expected. 
The shaded region on each graph shows the distribution 
of D for the minima directly connected to the global min- 
imum. Each of these minima is surrounded by its own 
basin of attraction, and a barrier must be surmounted be- 
fore the configuration point enters the catchment basin 
of the global minimum. Even so, in Fig. ||b for LJ13, 
the tail of the high-energy quench distribution overlaps 
somewhat with the distribution of connected minima, in- 
dicating that some of the points in the basin of attraction 
of the global minimum are as far from it as the closest 
local minima. In Fig. 0a for M13 with p = 4, however, 



the high-energy quench distribution overlaps completely 
with the distribution of connected minima. The overlap 
means that many configurations which differ structurally 
from the global minimum as much as the connected local 
minima still quench to the global minimum. Such con- 
figurations may include structures which for a slightly 
shorter-ranged potential fall into the catchment basin of 
a different local minimum. 
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FIG. 2. Distribution of Euclidean distances to the global 
minimum from configuration points in its basin of attraction 
for (a) M13, p = 4 and (b) LJ13. For each cluster, the distri- 
butions from microcanonical MD simulations at two energies 
are shown. At the lower energy, the majority of quenches 
lead to the global minimum, and at the higher energy, about 
half do so. The quench interval was 15 reduced time units. 
The duration of the high-energy simulations was 6 x 10 , and 
that of the low-energy ones was 3 x 10 s . The shaded areas 
show the distance distribution of minima that are directly 
connected to the global minimum. 

Further evidence to support this explanation comes 
from Fig. [I]. The MC results for p — 6 show that the in- 
crease in gradient of the caloric curve — indicative of the 
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system sampling a new region of configuration space — 
occurs at the temperature where the trajectory begins to 
escape from the global minimum, as shown by the de- 
crease in -Pgmi n - The transition is therefore out of the 
basin of attraction of the global minimum, effectively 
from solid-like to liquid- like states. The analogous re- 
sults for p = 4 show that the transition feature in the 
caloric curve starts before the probability of the global 
minimum drops significantly. This behaviour is sugges- 
tive of a weak transition from configurations close to the 
global minimum to higher-energy ones, all taking place 
within the basin of attraction of the global minimum it- 
self. 

The picture of the basin of attraction around the global 
minimum that emerges is therefore complex. It extends 



far into configuration space from the icosahedron, and 
past the catchment basins of local minima in certain 
directions. The structural dissimilarity of some points 
which arc formally associated with the global minimum 
suggests that quenching might not be the most mean- 
ingful way of dividing configuration space amongst the 
various minima. From this point of view, the underes- 
timation of the global minimum probability in Fig. |l|a 
by the harmonic superposition method is somewhat mis- 
leading, since it is not helpful to think of the distant con- 
figurations as belonging to the icosahedral well. The har- 
monic approximation may therefore give a more meaning- 
ful probability that the cluster has a structure resembling 
the global minimum than quenching. 
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FIG. 3. Relaxation of minima, grouped in "layers" away from the global minimum, for M13 at four values of the range 
parameter, p. In each plot, layer f is the global minimum, layer 2 contains all minima directly connected to layer I, etc. In 
each case, the microcanonical total energy is chosen such that the equilibrium probability of the global minimum is 0.4: for 
p = 4, 6, 10, 14, E/e = -33.17, -29.42, -28.43, -29.78, respectively. The inset for p = 14 shows the first half of the layer 5 
curve with a logarithmic time axis. The time is in units of (mrj/t) 1 ' 2 . 



In summary, although the harmonic approximation has 
only partial success in describing the equilibrium prop- 
erties of the clusters examined here, it is attractive in its 
simplicity, lack of empirical parameters, and clear physi- 
cal basis. More complicated analytic models do not nec- 



essarily provide greater insight, or even systematically 
improved results. We therefore adopt the harmonic ap- 
proximation for the rat e con stant and equilibrium prop- 
erty expressions in Sec. [IB, with the proviso that they 
will only be applied at low and moderate temperatures, 
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where the description should be adequate for our pur- 
poses. The resulting microcanonical rate constants, via 
Eq. (H), are given by 



kl(E) 



Fat pt(«-i) \E-Vi 



and the canonical expression from Eq. (|11[) is 



^PGt pt(« 



-(yt-VjO/fceT 



(21) 



(22) 



Finally, we note that the harmonic superposition method 
is likely to be more successful in the canonical ensemble 
than the microcanonical because at constant temperature 
the velocity distribution is independent of the configura- 
tion point. In contrast, at fixed total energy the kinetic 
energy is significantly further "above" the PES when the 
configuration point is near a deep potential well than 
when it is near higher-lying ones. 



B. Relaxation and the Range of the Potential 



assumption of stochastic transitions breaks down when 
the relaxation time approaches the vibrational period. 

The relaxation is straightforward for p — 4; the fur- 
thest layer decays whilst the populations of the global 
minimum and the intermediate layer grow monotonically. 
At p = 6 we see the accumulation and decay of a tran- 
sient population in layer 3 as the outward flow from these 
minima does not match the rapid initial input from the 
furthest layer. The layer 3 minima therefore constitute 
a kinetic bottleneck for relaxation down the funnel of 
the PES. More complex behaviour arises for p = 10 and 
14, where there are six layers of minima. The probabili- 
ties experience an initial jump as the system is released 
from its strongly non-equilibrium state, and then relax 
slowly to their final values. Layer 5 at p — 14 (shown 
in the inset of Fig. |J) shows particularly complicated be- 
haviour, rising suddenly at first, decaying slightly, and 
then rising again before relaxing monotonically. This os- 
cillation occurs because layer 4 develops a transient pop- 
ulation, blocking further downward output from layer 5 
while layer 6 is still releasing probability into layer 5 from 
above. 



1. Relaxation to the Global Minimum 



We now turn to the effect of the range of the poten- 
tial on the dynamics of the M13 clusters in the light of 
our previous analysis of the energy landscapes Q] . The 
larger number of rearrangements and smaller energy gra- 
dient on paths to the global minimum, as well as the 
higher downhill barriers, are expected to impede relax- 
ation to the global minimum as the range of the potential 
is decreased. 

Figure || illustrates the range dependence of struc- 
tural relaxation to the global minimum. The minima arc 
grouped into "layers" according to the smallest number 
of rearrangements required to reach the global minimum. 
Layer 1 contains only the global minimum, and level i + 1 
contains all minima directly connected to a minimum in 
level i, but not to any minimum in a layer lower than 
i. The initial probability vector at each value of p is 
a uniform distribution amongst the minima in the fur- 
thest layer from the global minimum. The energy has 
been chosen such that the equilibrium probability of the 
global minimum is 0.4, and the solution of the master 
equation is shown until the time at which 80% of this 
population (i.e., 0.32) has been achieved. 

The most striking trend is the increasing time scale as 
the potential range decreases, as predicted by the land- 
scape analysis [Q. The p = 14 cluster takes over three 
orders of magnitude longer to reach 80% of equilibrium 
than the p = 4 cluster. Note that for the latter system, 
the time scale plotted is only of the order of a few vibra- 
tional periods Q . In the light of previous work |ll[ the 
application of the master equation to the p — 4 cluster 
is therefore probably at the limit of validity, since the 



2. Search Times 



By analogy with the "folding time" in the protein fold- 
ing literature, we can examine the ability of the cluster 
to find its global minimum by defining a "search time" as 
the time taken for the probability of the global minimum 
to reach a particular value after the system is released 
from a non-equilibrium state. Here we use a probability 
threshold of 0.4, but we will discuss the effects of chang- 
ing this choice. 

Fig. H shows the search time for M13 as a function of 
temperature using four values of p. The initial probabil- 
ity vector was a uniform distribution amongst the minima 
in the layer furthest from the global minimum. The qual- 
itative shape of the curves is easily understood. As the 
temperature is increased from low values, the search time 
decreases because the thermal energy rises above the bar- 
riers between the minima. An optimal temperature, T opt 
is reached, where the search time is a minimum, r op t, 
above which it rises because the thermodynamic driving 
force towards the global minimum is reduced at higher 
temperatures. Ultimately, the equilibrium probability of 
the global minimum falls below the threshold of 0.4, and 
the search time is no longer defined. Similar behaviour 
has been observed for the same reasons in direct simula- 
tions of KC1 clusters |24|] and lattice protein models [g5| , 
as well as a master equation study of idealized energy 
landscapes 26 1. 
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FIG. 4. Search time as a function of temperature for M13 
at four values of the range parameter, p. The search time 
is defined as the time taken for the probability of the global 
minimum to reach 0.4, starting from an even probability dis- 
tribution amongst the minima in the layer furthest from the 
global minimum. The time is in units of (mrg/e) 1 / 2 . 



Table | lists the temperature and r op t , as well as the 
two temperatures, T\ ow and Thigh, at which the search 
time equals 2r op t. The difference Thigh — T\ ow provides 
a measure of the width of the temperature window for 
which searching is reasonably fast. The value of T opt in- 
creases with p, as expected from the relaxation profiles of 
the previous section, and this slowing down is accompa- 
nied by a decrease in the temperature width Thigh — T ow . 
As the range of the potential is decreased, the energy 
gap between the global minimum and the other minima 
becomes smaller, and the energy range spanned by the 
minima narrows. Hence, other minima come into play 
at lower temperature, and the temperature at which the 
global minimum ceases to dominate the equilibrium pop- 
ulations is lower when the potential is short-ranged. This 
observation explains why Thigh falls as p increases. For 
high p, the downhill barriers between minima are on av- 
erage larger, so that as the temperature is decreased, iso- 
merisation processes slow down more dramatically than 
for small p. Hence, the search time increases more rapidly 
as the temperature is lowered beyond T opt when the po- 
tential is short-ranged, resulting in the narrower ranges 
of Tiigh — T ow . 



TABLE I. "Searching" characteristics of M13 at four values of the range parameter, p. T opt is the optimal searching 
temperature, at which the search time is a minimum, T opt . Ti ow and Thigh (Ti ow < Thigh) are the two temperatures at which 
the search time equals 2r op t- 



p 


feTopt/e 




fcBTi ow /e 


ks Thigh /e 


Topt ^low 
Thigh Tiow 


4 


0.293 


1.47 


0.135 


0.344 


0.76 


6 


0.305 


6.54 


0.218 


0.333 


0.76 


10 


0.285 


130 


0.235 


0.301 


0.76 


14 


0.200 


2200 


0.170 


0.214 


0.68 



An analogy can be drawn here with the ease of fold- 
ing in proteins. The ratio, Tf/T g , of the "folding tem- 
perature" (where the native state becomes thermody- 
namically most stable) to the "glass transition temper- 
ature" (where the kinetics slow down dramatically) has 
been used as a measure of the ability of a protein to 
fold correctly p5| , p7t - Tf is roughly related to Thigh, 
which decreases with increasing p, whilst T g increases be- 
cause of the higher barriers for short-ranged potentials. 
Hence Tf/T g falls, in accordance with the observation 
that searching for the global minimum is harder when p 
is high. 

Interestingly, the curves in Fig. || for p = 4, 6 and 
10 have similar shapes. For example, the values of T op t 
differ by only 7%. Furthermore, T opt lies about three 
quarters of the way from T ow to Thigh in all three cases, 
as shown by the last column of Table |. p = 14 repre- 
sents an extreme case in which relaxation to the global 
minimum becomes very slow for values outside a small 
range near T opt . This optimum temperature is a com- 



promise between the slow dynamics at even moderately 
low temperatures, and the rapidly decreasing thermody- 
namic weight of the global minimum at moderately high 
temperatures. 

Choosing a different occupation probability of the 
global minimum as the criterion for the search time has 
predictable effects. If a higher threshold is used, the 
search time increases at any given temperature. The 
equilibrium probability of the global minimum drops be- 
low the threshold at a lower temperature, so the upper 
limit for which the search time is defined decreases. The 
search time rises more steeply below T opt because the 
probabilities must come closer to their equilibrium val- 
ues, and this approach is asymptotic. The combined ef- 
fect is that the search time curves all become narrower. 
However, for the cases tested (p — 4 and 6), T opt changed 
by only 10% as the threshold was varied from 0.2 to 0.5. 
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3. Relaxation of the Total Energy 

The evolution of the probability vector towards P cq 
is expressed macroscopically by the relaxation of some 
overall property, A, to its equilibrium value, A cq . If this 
property has a well defined value, Ai, for each state i in 
the master equation, the expectation value is a weighted 
average which can be expressed as a function of time us- 
ing Eq. (§): 



(A(t)) = AiPi® = £ 



c,e 



A,/ 



3=1 



3=2 



where 



i=i 



E^ } 

rn—l 



P m (0) 

Vp^ 



(23) 



(24) 



and the last expression in Eq. ( p3| ) uses the fact that 
the j = 1 term defines the baseline for relaxation, since 
Ai = 0. A mean relaxation time, r r , can be defined 
by normalizing the profile of (A) against t (such that it 
decays from 1 to 0), and evaluating the area under the 
resulting curve Q . For pure Debye (single exponential) 
relaxation, i.e., exp(— \t), one simply obtains r r = A -1 . 
Subtracting the equilibrium value from the right-hand 
side of Eq. (p3|), integrating from t = to oo, and nor- 
malizing using the value at t = yields 



E ^7 E 



(25) 



If the eigenvectors and eigenvalues of the transition ma- 
trix are not available, r r can be obtained by propagat- 
ing the master equation numerically until (A) has ef- 
fectively equilibrated, and then numerically integrating 
the normalized relaxation profile. Here we will examine 
the relaxation of the total energy of M13 as the popula- 
tions of the minima equilibrate at constant temperature. 
Within the harmonic approximation, therefore, we need 
Ai = Vi + nk B T in Eqs. @ and @. 

Theoretical treatments have shown that under par- 
ticular circumstances, the multi-exponential decay that 
arises from the master equation can lead to a variety 
of asymptotic behaviours. For example, the random en- 
ergy model, where the states have a Gaussian distribu- 
tion of energies, can lead either to stretched exponen- 
tial (oc exp[— (t/r) ]) or to power law (oc [t/£] p ) relax- 
ation of autocorrelation functions, depending on the form 
chosen for the transition rates |2£]. Palmer et al. also 
derived stretched exponential behaviour for a hierarchi- 
cally constrained model p(i| ]. In contrast Skorobogatiy 
et al. found power law and logarithmic (oc — hit) decay, 
but not stretched exponential decay, of the total energy 
in a protein model |pl| , depending on the temperature. 



Starting from a uniform distribution amongst the min- 
ima in the layer furthest from the global minimum, we 
found that none of the above forms (pure or stretched ex- 
ponential, power law, or logarithmic) gave a robust fit to 
the decay of the total energy from the master equation. 
At sufficiently long times, Eq. ( p3| ) approaches pure expo- 
nential behaviour, since all contributions except that of 
the slowest mode have decayed. At intermediate times, 
the parameters (particularly the stretching exponent, 9) 
that produced the best fit for the stretched exponential 
form were highly sensitive to the time interval over which 
data were supplied, and the resulting curves often devi- 
ated significantly from the master equation solution. The 
difficulty of obtaining an acceptable fit increased with de- 
creasing temperature, where the spread of the exponents 
Xj is wider. Of course, there is no reason why a general 
multi-exponential form like Eq. (^3|) should conform to 
any simplified model. 
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FIG. 5. Canonical ensemble Arrhenius plots for the re- 
laxation time of the total energy of M13 at (a) p = 4 and (b) 
p — 14. Circles are mean relaxation times from the master 
equation, dashed lines are fits to the Arrhenius form, and 
the solid line in (a) is a fit to the Vogel-Tammann-Fulcher 
(VTF) form. The relaxation time is in units of (mr c 2 / 'e) 1 ^ 2 . 
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Although the relaxation profiles are complicated, the 
mean relaxation times (integrated profil es) wer e found to 
follow simple empirical expressions. Fig. Ill B 3] shows the 
logarithm of the relaxation time as a function of inverse 
temperatur e for M 13 with p — 4 and 14. The p = 14 
plot in Fig. Ill B 3| b is well fitted by the Arrhenius form 
T r = r exp(A/fc B T) with T = 9.34 x 1CT 3 (mrg/e) 1 / 2 
and A = 2.00 e. The p = 4 plot in Fig. [ill B 3| a, however, 
shows significant deviation from the linearized Arrhenius 
expression. It is better fitted by the ubiquitous Vogel- 
Tammann-Fulcher (VTF) form [§2| 



T exp 



.4 



k B (T-T a ) 



(26) 



as shown by the solid line, for which the parameters are 
r = IMimrl/e) 1 /' 2 , A = 0.0699 e, k B T = 0.051 e. 
These values were found by least squares fitting of the 
logarithm of Eq. ( ^6|) with equally weighted points. 

The slower relaxation and larger database for the p = 
14 cluster mea nt that it was not computationally feasible 
to extend Fig. [IIB3b to lower temperature — the lowest 
point shown is at k B T = 0.13 e — and it is possible that 
deviation from the Arrhenius behaviour would occur be- 
low this value. However, the curvature of the p — 4 plot is 
clear over the same range, indicating that the relaxation 
dynamics of the two clusters respond differently to tem- 
perature changes. The origins of the difference probably 
stem from the decreasing slope of the energy landscape 
as the range of the potential is shortened. The energy 
intervals spanned by the minimum and transition state 
samples at higher p are narrower, making the landscape 
more uniform. This uniformity means that a change in 
the temperature has a similar effect on most of the indi- 
vidual interwell processes, each of which separately has 
an Arrhenius temperature dependence in the model we 
have used [Eq. (p2[)]. On the steeper landscape of the 
p = 4 cluster, however, there is a greater spread of lo- 
cal minimum energies and barrier heights, so that as the 
temperature is lowered, some processes become "frozen 
out" before others, resulting in longer relaxation times 
than expected by extrapolation of the high-temperature 
behaviour. 

In structural glasses, Arrhenius temperature depen- 
dence of relaxation times is associated with strong liq- 
uids, whereas VTF behaviour is indicative of fragility 
|j33f . If this classification can be applied to clusters, the 
results of this section suggest that increasing the range of 
the potential introduces a degree of fragility. Stillinger's 
picture of stron g li quids having a "uniformly rough" en- 
ergy landscape |2q| is in line with our previous analysis 
j| , which showed that decreasing the range of the poten- 
tial lowers the overall slope of the PES, leaving a flat but 
rough landscape. 



4- Relaxation Modes 

Is the time scale of relaxation mostly determined by 
the slowest relaxation mode of the master equation, i.e., 
the least negative non-zero eigenvalue of the transition 
matrix, or is the process it describes relatively unimpor- 
tant? Eqs. (|2^) and (^4j) show that the contribution of 
a given mode to the relaxation of a global property de- 
pends both on the property and on the initial probability 
distribution. However, we can still probe the nature of 
the probability flow described by a particular relaxation 
mode by comparing the size and sign of the components 
in the corresponding transition matrix eigenvector. Con- 
sider Eq. (j^) for a particular value of j; mode j makes 
an important contribution to the probability evolution of 
minimum i if uf^ — or equivalently -J P° q u^ — is large in 
magnitude. The mode represents an overall flow between 
minima i and k if and have opposite signs. 

The extreme eigenvalues and eigenvectors of the transi- 
tion matrix can be obtained efficiently for large matrices 
using Lanczos iteration |3^| . Inspection of the compo- 
nents of the eigenvectors for M43 at different values of p 
and different temperatures reveals some general trends. 
The extreme modes (i.e., the slowest and fastest) describe 
probability flow between a small number of minima, typ- 
ically fewer than five. The fastest modes tend to be be- 
tween minima that are directly connected by transition 
states. In contrast, the slowest modes are between un- 
connected minima, and probability flow involves interme- 
diate minima. The slow modes tend to involve one highly 
populated minimum. Hence, if the initial probabilities of 
the other minima that feature in these relaxation modes 
(with eigenvector components of opposite sign) are far 
from their equilibrium values, the slow modes may limit 
the overall relaxation. 

The number of minima participating in relaxation 
mode j can be measured using the index 



(27) 



which varies from 1 to n m ; n . Figure ^ shows hj as a func- 
tion of the eigenvalue, Xj, for all the relaxation modes of 
M13 with p = 6 at two temperatures. As observed above, 
the number of minima involved in the fastest and slowest 
modes is small. Many intermediate modes only involve 
a small number of minima too, but the modes that de- 
scribe more global flow are all clustered in the centre of 
the (logarithmic) scale. This pattern is more pronounced 
at the higher temperature. 

In a typical application of the master equation, there- 
fore, the initial processes involve rapid equilibration be- 
tween small groups of minima that are adjacent in config- 
uration space. This is followed by wider probability flow 
between larger groups of minima, and finally slow ad- 
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justment of the population of a few minima via processes 
involving multiple rearrangements. 
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FIG. 6. Number of minima participating in a relaxation 
mode of the master equation as a function of the eigenvalue 
of the mode for M13, p = 6, at (a) feT = 0.15 e and (b) 
feT = 0.40 e. The mode with zero eigenvalue has been ex- 
cluded. A is in units of (e/mr 2 ) 1 / 2 . 



IV. INTERFUNNEL DYNAMICS IN LJ 38 

We now turn to the double-funnel energy landscape of 
the 38-atom Lennard- Jones cluster, LJ38. The potential 
energy is given by 35 
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(28) 



The LJ 38 cluster is too large for a complete catalogue 
of minima and transition states to be obtained. How- 
ever, we have previously performed a thorough charac- 
terization of the low-energy regions of the PES using 
a database of 6000 minima and 8633 transition states 
[pjpql. The resulting disconnectivity graph |5 36 clearly 
showed two separated regions of configuration space, each 
having the character of a funnel. The global minimum 
lies at the bottom of a small funnel associated with 
28 minima that are characterized by face-centred cubic 
packing. The larger funnel of 446 minima leads to the 
lowest-energy icosahedrally packed minimum. Since this 
secondary funnel accounts for a large volume of config- 
uration space, and because the liquid-like minima are 
structurally more similar to the icosahcdral minima than 
the fee ones, the icosahedral funnel is expected to act as 
a kinetic trap for relaxation from high-energy states to 
the global minimum. 
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where a is the pair separation at which Vij = 0, and e is 
the pair well depth. We will use a and e as the units of 
the quantities they measure, setting both equal to unity; 
the topology of the PES is not affected by the values of 
these parameters. 



FIG. 7. Low-temperature properties of LJ38 calculated 
using the harmonic superposition approximation, (a) Equi- 
librium occupation probability of the fee and icosahedral fun- 
nels and the rest of configuration space, and (b) the heat 
capacity. 
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By grouping together the minima in each of the two 
funnels, we can study not only equilibration within the 
funnels as for M13, but also the dynamics between the 
funnels. Hence, we define the probabilities 

Pioc(t) = P *(*) aild P icos{t)= P *(*)' ( 29 ) 



The database obtained in previous work ||[36[ is a good 
representation of the low-lying regions of the PES, but 
does not extend far into the liquid-like regime, and so we 
are restricted to studying the dynamics at low tempera- 
tures, where the role of the liquid is less important. This 
is not a serious restriction, since the time scale separation 
of inter- and intrafunnel processes should be greatest at 
low temperature. 

To obtain an impression of the temperature range over 
which valid conclusions can be drawn, Fig. |?] presents 
some properties calculated using the harmonic superpo- 
sition approximation and the full low-energy database of 
6000 minima and 8633 transition states. Fig. ^a shows 
P^. and P^ s as a function of temperature. The global 
minimum only dominates at very low temperatures, since 
the larger number of minima and lower vibrational fre- 
quencies in the icosahedral funnel cause the phase volume 
of the latter to rise rapidly. At sufficiently high temper- 
ature the cluster should melt, and the occupation proba- 
bilities of the two funnels, which contain predominantly 
solid-like structures, should approach zero. However, the 
steady rise of the curve marked "rest" (i.e., 1— -Pf cc — Picos) 
is interrupted at k-&T ps 0.2 e, reflecting the fact that the 
sample of minima does not represent the relevant regions 
of the PES very well above this temperature. Figure uh 
shows the heat capacity, derived from Cy = (dU/dT)y 
and Eq. ([Tgj) . The small peak at k^T k, 0.12 e results 
from the transition from the fee to icosahedral regions of 
configuration space, and the main peak at k^T ss 0.18 e 
signifies the melting transition. These features are largely 
in agreement with a more sophisticated anharmonic su- 
perposition method, which was designed to model the 
thermodynamics from a representative sample of min- 
ima |2^j3^]. The harmonic results presented here pre- 
dict a lower peak for the melting transition, and beyond 
the melting temperature the heat capacity returns to its 
solid-like value. This result again reveals the deficiencies 
of the sample of minima in the liquid-like regime. We 
conclude, however, that as far as thermodynamic proper- 
ties are concerned, the present sample should be adequate 
for kBT < 0.2 e. The corresponding limit in the micro- 
canonical ensemble is E < —150 e, and the fee and icosa- 
hedral funnels have equal probability at E = — 160.5 e. 
We note that conventional simulations would not be able 
to measure the quantities in Fig. [7 reliably at such tem- 
peratures because the interfunnel dynamics are too slow. 
One of the aims of this section is to quantify the rate of 
passage between the funnels. 




\n[-X] 

FIG. 8. Spectra of transition matrix eigenvalues for LJ38 
in the microcanonical ensemble at (a) E — — 160 e and (b) 
E = —150 e. In each case, the eigenvalue of the most promi- 
nent interfunnel mode is marked by an arrow, showing the 
greater separation of its time scale from the other relaxation 
modes at low energy. A is in units of (e/mcr 2 ) 1 ^ 2 . 



A. Pruning the Database 

Because the time scale of interfunnel processes is so 
long at low temperatures, numerical integration of the 
master equation is not feasible, so we must diagonalize 
the rate matrix and use Eq. (|5|). However, diagonalisa- 
tion routines are likely to run into numerical problems 
when the matrix elements span many orders of magni- 
tude. Fig. H shows that lowering the temperature widens 
the spread of eigenvalues. At sufficiently low temper- 
ature, a number of eigenvalues begin to appear to be 
slightly positive, or the diagonalisation routines may fail 
altogether. Czerminski and Elber reported similar prob- 
lems H , and therefore restricted their studies to suffi- 
ciently high temperatures. In Sec. Ill B| we saw that the 
eigenvalues at the extremes of the spectrum tend to be 
associated with probability flow between a small number 
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of minima. If these minima are kinetically isolated, they 
cause the transition matrix to become nearly decompos- 
able, giving rise to the numerical difficulties. Since they 
do not participate in probability flow we should seek to 
exclude them. If they are not kinetically isolated, they 
will also appear in other, faster, relaxation modes. 

From a practical point of view, even if the transition 
matrix can be diagonalized without numerical difficulty, 
it is worth considering whether all the minima play a 
significant role. If some can be discarded without affect- 
ing the relaxation, the size of the matrix can be reduced, 
saving computational effort. Because of the astronomical 
number of minima on the PES, a large proportion of the 
minima found in a partial search are only linked to one 
other minimum in the database. Within the restricted 
sample, these minima constitute "dead ends" for proba- 
bility flow. They can act as buffers, absorbing and releas- 
ing probability as it flows towards equilibrium through 
the connected minimum, but they cannot act as path- 
ways for flow between minima or larger regions of the 
PES. The dead-end minima in our sample tend to be 
high in energy since the search algorithm only explores 
connections from low-energy minima thoroughly. When 
a high-energy minimum is found, the upward move is 
rejected and the search continues from the original mini- 
mum, never returning to the high ones. When modelling 
the probability flow from the bottom of one funnel into 
the other, high-energy dead-end minima are unlikely to 
play an important role because their equilibrium proba- 
bility is low and they do not mediate interfunnel flow. 

These considerations suggest ways of pruning the 
database. First, dead-end minima were identified, and 
were found to constitute about 70% of the sample of 6000 
minima. To eliminate the kinetically isolated minima at 
a given energy, the total outward rate constant was cal- 
culated for each dead-end minimum, and the minimum 
was discarded if the rate fell below a certain threshold. 
For example, at a total energy of E = — 160 e, the rate 
constants for individual processes span the range 10~ 41 
to 10° (e/mcr 2 ) 1 / 2 . Choosing a threshold of 10 12 , which 
corresponds to a time scale of seconds for argon parame- 
ters, reduces the sample of minima to 5944. A threshold 
of 10 -10 reduces it to 5861. While this may be suffi- 
cient to remove numerical difficulties in the diagonalisa- 
tion procedure, the matrices are still rather large. 

Trimming the sample on the basis of equilibrium prob- 
abilities reduced the number of minima more dramati- 
cally. At E — —160 e, discarding dead-end minima whose 
equilibrium probability was less than 10 -10 left just 1825 
minima, and a threshold of 10 -8 left 1782. Clearly, the 
number of minima removed by such a method decreases 
with increasing energy, since higher-energy states then 
become more populated. We will gauge the effect of 
pruning the database by examining the sensitivity of the 
results to the choice of the threshold. 



B. Interfunnel Rate Constants 

What is the rate of crossing between the fee and icosa- 
hedral funnels? We have previously shown that intercon- 
version of the fee and icosahedral minima is a multiple- 
step process [EJ — the lowest-energy path in our sample 
involves 13 successive rearrangements — but let us con- 
sider the overall scheme 



fee 



icos, 



(30) 



with "forward" and "reverse" rate constants k + and k_. 
The rate of change of the occupation probability of the 
fee funnel is accordingly 



dP icc (t) 
dt 



= -k+P {cc {t) + k-P icos (t). 



(31) 



From Fig. [?], we see that at sufficiently low temperature 
only the two funnels are significantly occupied at equilib- 
rium, as opposed to minima associated with the liquid- 
like state. Assuming that this is also the case away from 
equilibrium, provided that the initial probability is itself 
confined to the funnels, we can write 



PfccW+iW*) = 1. 
Using Eq. (J32|) and the equilibrium relationship 



k , P eq 



k P eq 



(32) 



(33) 



integration of Eq. (31) gives the basic result of first order 
kinetics for a two-state model: 



In 



Pfcc(t) 



-'fee 



Pfcc(0) 



poq 
^fec 



-(k++k-)t. 



(34) 



Figure ||a shows plots of Eq. and the analogous 
expression for Pi cos (t) in the microcanonical ensemble at 
E = — 160 e, starting from the global minimum. The 
plots were obtained from the analytic solution of the mas- 
ter equation after removing all dead-end minima with an 
equilibr ium p robability of less than 10 -8 , as described 



in Sec. IV A. The two lines are straight and coincide, 
and the slope yields k + + k- = 4.99 x 10~ 12 (e/mcr 2 ) 1 / 2 . 
This value closely matches one of the eigenvalues of the 
transition matrix, |A 4 | = 4.98 x 10~ 12 (e/mcr 2 ) 1 / 2 , sug- 
gesting that the corresponding eigenvector describes flow 
between the two funnels. The "net flow index" into or 
out of a funnel F produced by relaxation mode i of the 
master equation can be obtained by summing the com- 
ponents of eigenvector i that correspond to the minima 
belonging to F §§: 



ff = 



(35) 



If mode i represents probability flow between the funnels, 
then // cc and fl° os will be larger in magnitude than for 
other eigenvectors, and will have opposite signs, so that 
an increasing contribution is made to one funnel and a 
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decreasing contribution to the other, depending on the 
initial probability vector [see Eq. (§)]. At E = — 160 e, 
for the mode whose eigenvalue matches — (k + + we 
find / 4 fcc = 0.495 and f 4 cos = -0.495. The next largest 
net flow index for the fee funnel was /g| c = 3.97 x 10~ 4 
(with a corresponding index for the icosahedral funnel of 
f£° s = 9.94 x 10~ 7 ) and that for the icosahedral fun- 
nel was /g cos = —0.143 (with a corresponding index for 
the fee funnel of /| cc = 4.38 x 10~ 5 ). This result un- 
ambiguously identifies the fourth mode with interfunncl 
relaxation at this energy. 




-6 1 ■ ' ■ ' ■ 1 

10000 20000 30000 

time 

FIG. 9. Plots of Eq. (ji^) for LJ38 in the microcanonical 
ensemble with the initial probability of the global minimum 
set to unity. Solid lines are Eq. (B4|) for the fee funnel, and 
dashed lines are the equivalent for the icosahedral funnel, (a) 
E = -160 e, (b) E = -150 e. In (a) the lines coincide. The 
units of time are (a) 10 12 (ma 2 /e) 1/2 , and (b) (ma 2 /e) 1/2 . 

To test the effect of having pruned the database, we 
lowered the threshold for removal of minima from an 
equilibrium probability of 10~ 8 to 10~ 12 , resulting in a 
larger remaining sample of 2289 minima. The net flow 
index picked out relaxation mode 17, whose eigenvalue 
was A17 = —4.98 x 10~ 12 (e/mcr 2 ) 1 / 2 , i.e., the same as 
was obtained with the higher threshold. Replotting Fig. 



also yielded the same result as before. We note that, 
even though this new sample is less than half the size 
of the full database, the lowest eigenvalues already clash 
with the precision of the "zero" eigenvalue. Hence, some 
pruning is essential for numerical tractability and desir- 
able for computational speed, and it does not affect the 
result in this case. 

At low energy, the interfunnel mode is easily identified. 
As the energy is raised and more processes become "un- 
frozen" , the distinction is somewhat less clear. Figure 
^b shows Eq. ( |34| ) starting from the global minimum at 
E = — 150 e. The database was pruned using an equilib- 
rium probability threshold of 10 -8 and contained 3789 
minima. Whilst the decay of Pf cc (t) obeys the linearized 
relationship, the evolution of Pi cos (i) deviates from it in- 
creasingly as time progresses. This deviation is partly be- 
cause minima outside the two funnels have non-negligible 
populations, so Eq. ( |32"| ) does not hold, and also because 
-Picos(i) does not rise monotonically to its equilibrium 
value, but overshoots slightly and then decays. The net 
flow index picks out A95 = —9.19 x 10~ 5 (e/mcj 2 ) 1 / 2 with 
flf = 0.0859 and /^ c 5 os = -0.0686. These values are 
considerably smaller than those obtained at E = — 150 e. 
Although other modes may have a higher flow index for 
one funnel, the value for the other funnel is then either 
much smaller (indicating that the mode describes flow 
between one funnel and the non- funnel states), or of the 
same sign (indicating that flow is not between the fun- 
nels) . 

The slope of the solid line in Fig. ||b is —9.35 x 
10~ 5 (e/mc 2 ) 1 / 2 , which is not far from Ag5, but is actu- 
ally closer to A 96 = -9.39 x 10 ~ 5 (e/mcr 2 ) 1 / 2 . However, 
mode 96 is only weakly interfunnel: /|g c = 1.91 x 10~ 4 
and /gg os = —5.36 x 10~ 4 . At relatively high energies, 
where minima outside the funnels come into play and the 
simplified scheme of Eq. ( |30| ) breaks down, the net flow 
index therefore still provides a convenient way of iden- 
tifying the most important interfunnel relaxation mode 
and extracting the quantity (fc+ + k_). This method also 
has the advantage that it does not require evaluation of 
the master equation solution itself. 

The equilibrium relationship Eq. (|33|) allows the sepa- 
rate rate constants k + and fc_ for interfunnel flow to be 
obtained from the eigenvalue. Table || shows k + and fc_ 
as a function of temperature in the canonical ensemble. 
At each temperature, the database was pruned using an 
equilibrium probability threshold of 10 -8 , and the table 
shows how many of the full sample of 6000 minima and 
8633 transition states remained. At the lowest two tem- 
peratures, the pruning procedure removed all dead-end 
minima. The table also shows that the net flow indices 
for the interfunnel mode decrease in magnitude at the 
higher temperatures. This reduction is a result of the 
increasing involvement of higher-energy and non-funnel 
minima. Over the course of doubling the temperature, /ex- 
changes by ten orders of magnitude. For argon parame- 
ters, the span of time scales is hundreds of nanoseconds 
to an hour. 
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TABLE II. Interfunnel rate constants for LJ38 in the canonical ensemble. n' mln and n' ts are the numbers of minima 
and transition states remaining in the database after discarding minima with an equilibrium probability of less than 10~ 8 . 
/ fcc , J lcos and A are the net flow indices and eigenvalue of the interfunnel mode. A, k+, and fc_ are tabulated in units of 
(e/ma 2 ) 1/2 . 
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Figure |lOj shows that, over the temperature range in 
Table |n], k + and fc_ obey an Arrhenius temperature de- 
pendence law. Only very slight curvature is visible in the 
fe_ results. Fitting to the form k — Aexp(~E a /k-QT) 
gives the pre-exponential factors and activation energies 
for the forward and reverse processes: 

fc+: fcc ^ icos A = 11.1 (e/mer 2 ) 1/2 £ a = 4.12e 
icos ^ fcc A = 3.18 (e/mer 2 ) 1/2 £ a = 3.19e 
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FIG. 10. Arrhenius plot of k+ (circles) and fc_ (squares) 
for the interfunnel dynamics of L J38 . The lines are fits to the 
form k = Aexp(— E^/ksT). The units of the rate constant 
are (e/ma 2 ) 1/2 . 

Interestingly, the effective activation energy for fcc — > 
icos is close to the overall barrier on the lowest-energy 
path between the funnels, starting from the global min- 
imum ||, which is 4.22 e. This result suggests that the 
pathways passing over the highest-energy transition state 
on the lowest-energy pathway determine the interfunnel 
dynamics. E a for icos — > fcc, however, is significantly 
lower than the overall potential barrier for the reverse 



process, which is 3.54 e, starting from the lowest-energy 
icosahedral minimum. The discrepancy can be attributed 
to the fact that more than one minimum in the icosahe- 
dral funnel is substantially occupied. Hence, the effective 
barrier from the icosahedral funnel should not be mea- 
sured from the lowest-energy minimum in the icosahedral 
funnel, but with respect to a weighted average of the min- 
imum energies, X^ieicos PNi- The slight curvature in the 
Arrhenius plot for fc„ is a result of the temperature de- 
pendence of this average. 

An alternative method for computing the rates for in- 
terfunnel conversion would be to calculate first the free 
energy barrier between the two funnels, and then the 
transmission coefficient for passage over this barrier |}9| . 
The free energy barriers have been calculated for LJ38 
and were found to decrease as the temperature is in- 
creased towards that required for melting. Therefore, if 
this method is also to show an Arrhenius behaviour, the 
temperature dependence of the transmission coefficient 
must compensate that of the free energy barriers. 

We have seen that the interfunnel rates drop dramat- 
ically as the temperature is lowered. At the same time, 
the increasing net flow indices show that the correspond- 
ing relaxation modes become more distinct from other 
processes. These features concur with Stillinger's inter- 
pretation of a processes in fragile liquids ^8|. In this 
picture, (3 processes are faster and more localized in con- 
figuration space, whereas the a processes, which become 
relatively slow at low temperatures, are thought to oc- 
cur between "craters" (a similar concept to funnels) on 
the energy landscape. These differences give rise to a 
bifurcation of time scales, which is visible for LJ38 in 
the eigenvalue spectra of Fig. ||. However, a processes 
have been observed to have a non- Arrhenius tempera- 
ture dependence, in contrast with the results of Fig. [l(]. 
The Arrhenius behaviour may be a genuine feature of 
the dynamics of the LJ38 cluster, but could also be a 
limitation of the present approach. In particular, our 
incomplete database contains only a limited number of 
pathways between the two funnels, and competition be- 
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tween such paths would give rise to non-Arrhenius be- 
haviour. We can study only a limited temperature range 
with our database, and the modelling of individual pro- 
cesses might not be sufficiently sophisticated. We note 
that Angelani et al. also observed unexpected Arrhcnius 
behaviour in their master equation study of a fragile glass 
former [EoLEil. 



C. Equilibration 

The progress of the probability vector towards equi- 
librium can be visualized using equilibration graphs 
p6| , ^2| -|44[ . Such a graph has a time axis, on which lines 
denote a group of states in local equilibrium with each 
other. Nodes join lines at the time that the correspond- 
ing groups first come into equilibrium, until there is just 
one group and the whole system has equilibrated. We de- 
fine the time that minima i and j come into equilibrium 
as the smallest value of t after which 



cq I 



(36) 



is always satisfied, and in the present work we set £ = 
0.01. 

Figure [H] shows equilibration graphs for the six min- 
ima in each of the two funnels of LJ38 that have the 
greatest equilibrium probability at E = — 150 e. They 
are numbered in order of increasing probability within 
each funnel. The initial vertical position of each mini- 
mum is taken as the integrated path length, 5 gm , along 
the shortest path to the global minimum, making the 
two groups clearly distinguishable on this axis. Three 
microcanonical energies are plotted, spanning the range 
of applicability of our database. The sample was pruned 
using an equilibrium threshold of 10 -8 at each energy. 
The evolution of two initial probability vectors was con- 
sidered: in the left-hand graphs, the initial probability 
of the global minimum is unity, and in the right-hand 
graphs the probability commences in the lowest-energy 
icosahedral minimum. 

In each of the six graphs, the minima within a funnel 
come into equilibrium with each other before the sepa- 
rate funnels do so. This result explicitly demonstrates 
the longer time scale of the interfunnel dynamics. The 
order of equilibration within each funnel is the same at all 
three energies studied, irrespective of the funnel in which 
the probability is initiated. There is a small exception in 
graph (e), where minimum 6' in the icosahedral funnel 
equilibrates with 3' and 5' before it joins 1', 2' and 4'. 
Interestingly, the lowest-energy icosahedral minimum, 3', 
is one of the last minima to reach equilibrium within the 



icosahedral funnel, presumably because of the high bar- 
riers surrounding it (|] . 

The equilibration of the fee funnel is much more sensi- 
tive to the initial probability than that of the icosahedral 
funnel, in spite of the fact that the rate constants, fc+ 
and are roughly equal at E = — 160 e. This differ- 
ence in behaviour of the two funnels arises from the fact 
that the absolute probabilities of the minima in the fee 
funnel, other than the global minimum itself, are several 
orders of magnitude smaller than those of the minima in 
the icosahedral funnel. Hence, small changes in probabil- 
ity due to transient flows easily disturb the equilibrium 
between minima in the fee funnel, since distance from 
equilibrium is measured relative to the final probability 
by the left-hand side of inequality ([36]). When the cluster 
starts in the global minimum, the probability of the fee 
minima decreases monotonically, and the minima within 
the funnel can equilibrate with each other rapidly. When 
the probability is initialized in the lowest-energy icosahe- 
dral minimum, however, the influx of probability to the 
fee funnel must pass through the minima within the fun- 
nel on its way to the global minimum, and the transients 
must settle down almost completely before equilibrium is 
permanently established. This effect is reflected by the 
shift of the equilibration nodes of the fee funnel to later 
times as one goes from the left-hand to right-hand side 
equilibration graphs in Fig. 

At sufficiently low energy, the global potential energy 
minimum must be the most populated state at equilib- 
rium. However, we have seen that there is a kinetic bot- 
tleneck to entering its funnel. Hence, if the cluster is 
prepared in a liquid-like state, it is most likely to col- 
lapse into the icosahedral funnel, which is larger and 
structurally more similar to the liquid, even though it 
is not the equilibrium state. Over a sufficiently long 
time, the cluster must then convert to the global min- 
imum. Our master equation model shows this behaviour 
clearly. Starting from a uniform distribution amongst 
the 25 highest-energy minima in the sample, P{ cc (t) and 
Picos(t) were monitored as the system evolved towards 
equilibrium at the low energy of E — —160 e. The results 
are shown in Fig. O. 

The initial states decay rapidly into other non-funnel 
states, and the two funnels experience a slow increase 
in population. Although some probability enters the 
fee funnel, it reaches a plateau while the population of 
the icosahedral funnel continues to grow. This growth 
reaches a maximum before eventually decaying towards 
its equilibrium value, as the probability trickles into the 
global minimum. The icosahedral funnel acts as a kinetic 
trap, and only releases the cluster into the global mini- 
mum on a long time scale. We note that direct simulation 
of the trapping effect by standard MD would therefore be 
highly problematic. 
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FIG. 11. Equilibration graphs for LJ38 in the microcanonical ensemble at three energies: (a) and 
(b) — 160 e, (c) and (d) — 155 e, (e) and (f) — 150 e. In each row, the left-hand graph is for the global 
minimum (labelled 1) having an initial probability of 1, and the right-hand graph is for the lowest-energy 
icosahedral minimum (labelled 3') having an initial probability of 1. Lines representing individual minima 
commence at a vertical position corresponding to the shortest integrated path length, S gln , to the global 
minimum, and an arbitrary horizontal position. Nodes join lines when the corresponding states first come 
into equilibrium. In (a), unprimed numbers indicate minima in the fee funnel, and primed numbers indicate 
minima in the icosahedral funnel. The time is in units of (ma 2 /e) 1 ^ 2 . 
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FIG. 12. Relaxation of LJ38 from high-energy minima at 
a total energy of — 160 e, showing the fast and slow contri- 
butions to the final probability of the fee funnel. The time 
is in units of (ma 2 /e) 1 ^ 2 . 

Although Fig. [l2| unambiguously demonstrates the 
separate fast and slow contributions to the relaxation, 
the precise partitioning between the two funnels at the 
plateau stage can depend on the initial probability distri- 
bution. The distribution chosen here is rather artificial 
because our sample of minima does not extend into the 
liquid-like range. If we could release the system from a 
high-temperature liquid-like state, the fraction of proba- 
bility flowing into the icosahedral funnel would probably 
be even larger because of the greater structural similarity 
of the liquid-like structures with the icosahedral rather 
than fee minima. Although the fraction of fee minima is 
already small in our database, it would be much smaller 
in a more comprehensive sample. 

Two-stage dynamics have been observed experimen- 
tally in the folding of hen egg lysozyme, in which two 
routes to the native state have been postulated: one fast 
and direct, the other passing via partially folded con- 
formations which act as kinetic traps and reorganize to 
the native state only slowly p5[ . Conversely, the pro- 
tein plasminogen activator inhibitor 1 rapidly folds to 
the active state, but converts to an inactive form on a 
much slower time scale implying that the active 

state is not the global free energy minimum but is only 
metastable, like the icosahedral minima in LJ38 at low 
energy. Similarly, human prion protein has two long- 
lived forms, whose relative stability can be adjusted by 
varying the pH. The time scale for conversion to the more 
stable form is of the order of days p7[ . 



longer than those accessible by direct simulation, and de- 
scribes the behaviour of an equilibrating ensemble with- 
out the need to average over separate trajectories. 

The harmonic approximation for the density of states 
of individual minima and transition states provides a sim- 
ple but physically clear basis for calculating equilibrium 
properties and rate constants. Provided it is not applied 
at excessively high temperatures, it gives a qualitatively 
useful description of the thermodynamic and dynamic 
properties of the energy landscape. 

As we predicted Q] , relaxation to the global minimum 
is slower when the range of the potential is shorter. An 
optimal temperature for this relaxation is obtained by a 
compromise between the decreasing rates at low temper- 
atures and the decreasing thermodynamic driving force 
at high temperatures. When the range of the potential 
is long, the cluster exhibits a wide temperature window 
over which relaxation is quite efficient. In contrast, when 
the range is short, small deviations from the optimal tem- 
perature hinder the rate appreciably. 

Although the relaxation profiles of the total energy at 
fixed temperature do not appear to be well described by 
any of the commonly used empirical forms, the tempera- 
ture dependence of the mean relaxation time followed an 
Arrhenius law for p — 14 and a Vogel-Tammann-Fulchcr 
law at p — 4. These results again reflect the greater uni- 
formity of the short-range PES that was deduced in the 
landscape analysis. 

Application of the analytic solution of the master equa- 
tion to the low-energy database of LJ38 required the re- 
moval of unimportant minima and transition states from 
the sample. "Dead-end" minima were removed if their 
equilibrium probability fell below a low threshold at the 
temperature of interest. The relaxation modes of the 
resulting database were analysed using a flow index to 
extract the rate of passage between the two funnels on 
the energy landscape at low temperatures. The equili- 
bration patterns within and between the funnels clearly 
revealed the double-funnel structure. High-energy distri- 
butions relaxed preferentially into the secondary funnel 
of icosahedral minima rather than the close-packed fun- 
nel surrounding the global minimum. This behaviour 
stems from the greater structural similarity of the liquid 
to the icosahedral minima, which is reflected in the pat- 
terns of connectivity on the PES. Eventually, the cluster 
escaped from this kinetic trap into the global minimum, 
which is thermodynamically favoured at sufficiently low 
temperature. 



V. SUMMARY 
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